home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 49 / Amiga Format CD49 (2000-01-17)(Future Publishing)(GB)(Track 1 of 3)[!][issue 2000-02].iso / -serious- / programming / c / tkalgebra / lib / uklad.cpp < prev    next >
Text File  |  1999-12-06  |  10KB  |  590 lines
  1. #include "uklad.h"
  2. #include "macierz.h"
  3. #include <math.h>
  4. #include <stdio.h>
  5. #include <time.h>
  6.  
  7. double UkladNO::Dokladnosc=10E-100;
  8. char UkladNO::Czy_Zmieniac=1;
  9. double mymax(double a, double b)
  10. {  if(fabs(a)>fabs(b))
  11.         return fabs(a);
  12.     else
  13.         return fabs(b);
  14. }
  15.  
  16.  
  17. UkladSt::UkladSt(Macierz &t,Wektor &z)
  18. {  x=NULL;
  19.     pamiec=0;
  20.     clock_t sta,end;
  21.     A=&t;
  22.     ERR=0;
  23.     LU=new Macierz(t);
  24.     b=&z;
  25.     Czas=0;
  26.     Bledy=BEZ_BLEDU;
  27.     if(LU)
  28.     {  if(LU->Bledy)
  29.             Bledy=BLEDNE_DANE|LU->Bledy;
  30.         else
  31.         {  Bledy=BEZ_BLEDU;
  32.             sta=clock();
  33.             if(!(LU->GaussToLU()))
  34.                 Bledy|=BLEDNE_OPERACJE;
  35.             else
  36.             {  if(!Rozwiaz())
  37.                     Bledy=ZLE_UWARUNKOWANIE;
  38.             }
  39.             end=clock();
  40.             Czas=(double)(end-sta)/CLOCKS_PER_SEC;
  41.         }
  42.     }
  43.     else
  44.         Bledy=BRAK_PAMIECI;
  45. }
  46.  
  47.  
  48. UkladSt::~UkladSt()
  49. {  if(x)
  50.         delete x;
  51.     x=NULL;
  52.     if(LU)
  53.         delete LU;
  54.     LU=NULL;
  55.     if(pamiec&&A)
  56.         delete A;
  57.     if(pamiec&&b)
  58.         delete b;
  59. }
  60.  
  61. int UkladSt::Rozwiaz()
  62. {  if(LU->Bledy||Bledy)
  63.     {  Bledy|=LU->Bledy;
  64.         return 0;
  65.     }
  66.     else
  67.     {  x=Rozwiazanie(*LU,*b);
  68.         if(x)
  69.         {  if(x->Bledy)
  70.             {  Bledy|=BLEDNE_OPERACJE;
  71.                 return 0;
  72.             }
  73.             else
  74.             {  ERR=SprERR();
  75.                 return 1;
  76.             }
  77.         }
  78.         else
  79.         {  Bledy|=BRAK_PAMIECI;
  80.             return 0;
  81.         }
  82.     }
  83.  
  84. }
  85.  
  86. void UkladSt::Print(const char* nazwa)
  87. {  FILE *plik;
  88.     plik=fopen(nazwa,"w+");
  89.     if(plik)
  90.     {  fprintf(plik,"Macierz A:\n");
  91.         A->Print(plik);
  92.         fprintf(plik,"Wektor b:\n");
  93.         b->Print(plik);
  94.         fprintf(plik,"Rozwi±zanie:\n");
  95.         x->Print(plik);
  96.     }
  97.     fclose(plik);
  98. }
  99.  
  100.  
  101. UkladNO::UkladNO()
  102. {  Bledy=0;
  103.     ERR=0;
  104.     pamiec=0;
  105.     x=NULL;
  106.     A=NULL;
  107.     b=NULL;
  108.     I_K=NULL;
  109.     Q=NULL;
  110.     R=NULL;
  111.     y=NULL;
  112.     B=NULL;
  113.     d=NULL;
  114.     Czas=0;
  115. }
  116.  
  117. UkladNO::UkladNO(Macierz &t,Wektor &z)
  118. {  x=NULL;
  119.     clock_t sta,end;
  120.     Bledy=0;
  121.     ERR=0;
  122.     A=NULL;
  123.     b=NULL;
  124.     I_K=NULL;
  125.     Q=NULL;
  126.     R=NULL;
  127.     y=NULL;
  128.     B=NULL;
  129.     d=NULL;
  130.     Czas=0;
  131.     A=&t;
  132.     int i;
  133.     pamiec=0;
  134.     No0_K=A->Ile_Kolumn;
  135.     Bledy=BEZ_BLEDU;
  136.     b=&z;
  137.     Q=new Macierz(t);
  138.     R=new Macierz(Q->Ile_Kolumn,Q->Ile_Kolumn);
  139.     y=new Wektor(Q->Ile_Kolumn);
  140.     B=new Wektor(z);
  141.     I_K=new unsigned int[A->Ile_Kolumn];
  142.     if(B&&Q&&I_K)
  143.     {  //A->Print();
  144.         for(i=0;i<A->Ile_Kolumn;i++)
  145.             I_K[i]=i;
  146.         sta=clock();
  147.         if((t.Bledy==BEZ_BLEDU)&&(z.Bledy==BEZ_BLEDU))
  148.         {  if(!(OrtogonalizacjaGS()))
  149.                 Bledy|=ZLE_UWARUNKOWANIE;
  150.             else
  151.             {  ZnajdzB();
  152.                 if(!Rozwiaz())
  153.                     Bledy=ZLE_UWARUNKOWANIE;
  154.             }
  155.         }
  156.         else
  157.             Bledy|=BLEDNE_DANE;
  158.         end=clock();
  159.         Czas=(end-sta)/CLOCKS_PER_SEC;
  160.     }
  161.     else
  162.         Bledy=BRAK_PAMIECI;
  163. }
  164.  
  165. UkladNO::~UkladNO()
  166. {  if(Q)
  167.         delete Q;
  168.     if(R)
  169.         delete R;
  170.     if(y)
  171.         delete y;
  172.     if(I_K)
  173.         delete []I_K;
  174.     if(B)
  175.         delete B;
  176.     if(x)
  177.         delete x;
  178.     if(d)
  179.         delete d;
  180.     d=NULL;
  181.     x=NULL;
  182.     if(pamiec&&A)
  183.         delete A;
  184.     if(pamiec&&b)
  185.         delete b;
  186. }
  187.  
  188. int UkladNO::ZamienK(int ktora)
  189. {  unsigned int i;
  190.     int zk=ktora;
  191.     Wektor *PX;
  192.     double maxnorm,inf;
  193.     PX=(Q->Kolumny[I_K[ktora]]);
  194.     maxnorm=PX->Norma2();
  195.     inf=maxnorm;
  196.     double b;
  197.     for(i=ktora+1;i<No0_K;i++)
  198.     {  PX=(Q->Kolumny[I_K[i]]);
  199.         b=PX->Norma2();
  200.         if(maxnorm<b)
  201.         {  maxnorm=b;
  202.             zk=i;
  203.         }
  204.     }
  205.     if(inf==0)
  206.     {   No0_K--;
  207.     }
  208.  
  209.     if(zk>ktora)
  210.     {  i=I_K[ktora];
  211.         I_K[ktora]=I_K[zk];
  212.         I_K[zk]=i;
  213.         return 1;
  214.     }
  215.     else
  216.     {  if(maxnorm<Dokladnosc)
  217.         {  if(ktora>0)
  218.                 No0_K=ktora-1;
  219.             else
  220.                 No0_K=0;
  221.             return 0;
  222.         }
  223.         else
  224.             return 1;
  225.     }
  226. }
  227.  
  228.  
  229. int UkladNO::OrtogonalizacjaGS()
  230. {  int i,j,k;
  231.     d=new Wektor(Q->Ile_Wierszy+1);
  232.     if(R&&y)
  233.     {  for(i=0;i<R->Ile_Wierszy;i++)
  234.         {  for(j=0;j<i;j++)
  235.                 (*R)[i][j]=0;
  236.             (*R)[i][i]=1;
  237.         }
  238.         for(k=1;k</*A->Ile_Kolumn*/No0_K+1;k++)
  239.         {  if(Czy_Zmieniac==1)
  240.                 ZamienK(k-1);
  241.             (*d)[k]=0;
  242.             for(i=0;i<Q->Ile_Wierszy;i++)
  243.                 (*d)[k]+=((*Q)[i][I_K[k-1]]*(*Q)[i][I_K[k-1]]);
  244.             if((*d)[k]==0)
  245.             {  No0_K=i;
  246.                 return i;
  247.             }
  248.             for(j=k;j<No0_K;j++)
  249.             {  (*R)[k-1][j]=0;
  250.                 for(i=0;i<Q->Ile_Wierszy;i++)
  251.                     (*R)[k-1][j]+=((*Q)[i][I_K[k-1]]*(*Q)[i][I_K[j]]);
  252.                 (*R)[k-1][j]=((*R)[k-1][j]/(*d)[k]);
  253.                 for(i=0;i<Q->Ile_Wierszy;i++)
  254.                     (*Q)[i][I_K[j]]-=((*R)[k-1][j]*(*Q)[i][I_K[k-1]]);
  255.             }
  256.             
  257.         }
  258.     }
  259.     else
  260.         return 0;
  261.     return 1;
  262. }
  263.  
  264. void UkladNO::ZnajdzB()
  265. {  int i,k;
  266.     for(k=1;k<No0_K+1;k++)
  267.     {  (*y)[k-1]=0;
  268.         for(i=0;i<Q->Ile_Wierszy;i++)
  269.             (*y)[k-1]+=((*Q)[i][I_K[k-1]]*(*B)[i]);
  270.         (*y)[k-1]/=(*d)[k];
  271.         for(i=0;i<Q->Ile_Wierszy;i++)
  272.             (*B)[i]-=((*R)[0][k-1]*(*Q)[i][I_K[k-1]]);
  273.     }
  274. }
  275.  
  276. int UkladNO::Rozwiaz()
  277. {  int i,j;
  278.     double pomoc;
  279.     if(Q&&B)
  280.     {  if((Q->Bledy==BEZ_BLEDU)&&(B->Bledy==BEZ_BLEDU))
  281.         {  if((R->Bledy==BEZ_BLEDU)&&(y->Bledy==BEZ_BLEDU))
  282.             {  x=new Wektor(R->Ile_Wierszy);
  283.                 if(x)
  284.                 {  if(x->Bledy==BEZ_BLEDU)
  285.                     {  for(i=No0_K;i<R->Ile_Kolumn;i++)
  286.                             (*x)[I_K[i]]=0;
  287.                         for(i=No0_K-1;i+1;i--)
  288.                         {  pomoc=0;
  289.                             for(j=/*R->Ile_Kolumn*/No0_K-1;j>i;j--)
  290.                                 pomoc+=((*R)[i][j]*(*x)[I_K[j]]);
  291.                             (*x)[I_K[i]]=(*y)[i]-pomoc;
  292.                         }
  293.                     }
  294.                     else
  295.                         Bledy|=BRAK_PAMIECI;
  296.                 }
  297.                 else
  298.                     Bledy|=BRAK_PAMIECI;
  299.             }
  300.             else
  301.                 Bledy|=BLEDNE_OPERACJE;
  302.         }
  303.         else
  304.             Bledy|=BLEDNE_DANE;
  305.     }
  306.     else
  307.         Bledy|=BRAK_PAMIECI;
  308.     if(Bledy==BEZ_BLEDU)
  309.     {  ERR=SprERR();
  310.         return 1;
  311.     }
  312.     else
  313.         return 0;
  314. }
  315.  
  316. void UkladNO::Print(const char* nazwa)
  317. {  FILE *plik;
  318.     plik=fopen(nazwa,"w+");
  319.     if(plik)
  320.     {  fprintf(plik,"Macierz A:\n");
  321.         A->Print(plik);
  322.         fprintf(plik,"Wektor b:\n");
  323.         b->Print(plik);
  324.         fprintf(plik,"Rozwi±zanie:\n");
  325.         x->Print(plik);
  326.     }
  327.     fclose(plik);
  328. }
  329.  
  330. Uklad *RozwiazanieUkladu(Macierz &A, Wektor &b, char Flaga)
  331. {  Uklad *U;
  332.     Macierz *Z,*pom;
  333.     Wektor *BB;
  334.     int i;
  335.     if((A.Bledy==BEZ_BLEDU)&&(b.Bledy==BEZ_BLEDU))
  336.     {  if((A.Ile_Kolumn==A.Ile_Wierszy)&&(A.Ile_Wierszy==b.Ile_Wierszy))
  337.         {  if((Flaga==AUTODETECT_)||(Flaga==GAUSS_))
  338.                 return new UkladSt(A,b);
  339.             else
  340.                 return new UkladNO(A,b);
  341.         }
  342.         else
  343.         {  if((A.Ile_Kolumn<A.Ile_Wierszy)&&(A.Ile_Wierszy==b.Ile_Wierszy))
  344.             {  if((Flaga==AUTODETECT_)||(Flaga==ORTOGONALIZACJAGS_))
  345.                     return new UkladNO(A,b);
  346.                 else
  347.                 {  pom=A.Transponuj();
  348.                     if(pom)
  349.                     {  Z=new Macierz();
  350.                         BB=new Wektor();
  351.                         (*Z)=(*pom)*A;
  352.                         (*BB)=(*pom)*b;
  353.                         delete pom;
  354.                         if((BB->Bledy==Z->Bledy)&&(Z->Bledy==BEZ_BLEDU))
  355.                         {  U= new UkladSt(*Z,*BB);
  356.                             U->pamiec=1;
  357.                             return U;
  358.                         }
  359.                         else
  360.                             return NULL;
  361.                     }
  362.                     else
  363.                         return NULL;
  364.                 }
  365.             }
  366.             else
  367.             {  if((A.Ile_Kolumn>A.Ile_Wierszy)&&(A.Ile_Wierszy==b.Ile_Wierszy))
  368.                 {  if((Flaga==AUTODETECT_)||(Flaga==ORTOGONALIZACJAGS_))
  369.                     {
  370.                         return new UkladPO(A,b);
  371.                     }
  372.                     else
  373.                         return NULL;
  374.                 }
  375.                 else
  376.                     return NULL;
  377.             }
  378.         }
  379.     }
  380.     else
  381.     {  return NULL;
  382.     }
  383. }
  384.  
  385. double Uklad::SprERR()
  386. {  double pom,ret=0;
  387.     SrBlad=0;
  388.     unsigned int i,j;
  389.     if(x)
  390.     {  if(x->Bledy==BEZ_BLEDU)
  391.         {  for(i=0;i<A->Ile_Wierszy;i++)
  392.             {  pom=0;
  393.                 for(j=0;j<A->Ile_Kolumn;j++)
  394.                     pom+=((*A)[i][j])*((*x)[j]);
  395.                 ret=mymax((pom-((*b)[i])),ret);
  396.                 SrBlad+=fabs(pom-(*b)[i]);
  397.             }
  398.             if(SrBlad>0)
  399.                 SrBlad=SrBlad/A->Ile_Wierszy;
  400.         }
  401.         else
  402.             ret=257;
  403.     }
  404.     else
  405.         ret=257;
  406.     return ret;
  407. }
  408.  
  409. double DokladnoscMO(Macierz &A,Macierz &W)
  410. {  int i,j;
  411.     Macierz I;
  412.     double err=0;
  413.     int xz;
  414.     if((A.Bledy==BEZ_BLEDU)&&(W.Bledy==BEZ_BLEDU))
  415.     {  I=A*W;
  416.         if(I.Bledy==BEZ_BLEDU)
  417.         {  for(i=0;i<I.Ile_Wierszy;i++)
  418.             {  for(j=0;j<I.Ile_Kolumn;j++)
  419.                 {  if(i==j)
  420.                         err=mymax(err,((I[i][j])-1));
  421.                     else
  422.                         err=mymax(err,I[i][j]);
  423.                 }
  424.             }
  425.         
  426.         }
  427.         else
  428.             err=256;
  429.     }
  430.     else
  431.         err=257;
  432.     return err;
  433. }
  434.  
  435. UkladPO::UkladPO(Macierz &t,Wektor &z)
  436. {  x=NULL;
  437.     Bledy=0;
  438.     clock_t sta,end;
  439.     ERR=0;
  440.     A=NULL;
  441.     b=NULL;
  442.     I_K=NULL;
  443.     Q=NULL;
  444.     R=NULL;
  445.     y=NULL;
  446.     B=NULL;
  447.     d=NULL;
  448.     RT1=NULL;
  449.     Czas=0;
  450.     A=&t;
  451.     int i;
  452.     char tmp;
  453.     pamiec=0;
  454.     Bledy=BEZ_BLEDU;
  455.     b=&z;
  456.     Q=t.Transponuj();
  457.     No0_K=Q->Ile_Kolumn;
  458.     R=new Macierz(Q->Ile_Kolumn,Q->Ile_Kolumn);
  459.     y=new Wektor(Q->Ile_Wierszy);
  460.     B=new Wektor(Q->Ile_Kolumn);
  461.     I_K=new unsigned int[Q->Ile_Kolumn];
  462.     if(B&&Q&&I_K)
  463.     {  for(i=0;i<Q->Ile_Kolumn;i++)
  464.             I_K[i]=i;
  465.         sta=clock();
  466.         if((t.Bledy==BEZ_BLEDU)&&(z.Bledy==BEZ_BLEDU))
  467.         {  tmp=Czy_Zmieniac;
  468.             Czy_Zmieniac=0;
  469.             if(!(OrtogonalizacjaGS()))
  470.                 Bledy|=ZLE_UWARUNKOWANIE;
  471.             else
  472.             {  ZnajdzB();
  473.                 if(!Rozwiaz())
  474.                     Bledy=ZLE_UWARUNKOWANIE;
  475.             }
  476.             Czy_Zmieniac=tmp;
  477.         }
  478.         else
  479.             Bledy|=BLEDNE_DANE;
  480.         end=clock();
  481.         Czas=(end-sta)/CLOCKS_PER_SEC;
  482.     }
  483.     else
  484.         Bledy=BRAK_PAMIECI;
  485. }
  486.  
  487. UkladPO::~UkladPO()
  488. {  if(Q)
  489.         delete Q;
  490.     Q=NULL;
  491.     if(R)
  492.         delete R;
  493.     R=NULL;
  494.  
  495.     if(y)
  496.         delete y;
  497.     y=NULL;
  498.     if(I_K)
  499.         delete []I_K;
  500.     I_K=NULL;
  501.     if(B)
  502.         delete B;
  503.     B=NULL;
  504.     if(x)
  505.         delete x;
  506.     if(RT1)
  507.         delete RT1;
  508.     x=NULL;
  509.     if(pamiec&&A)
  510.         delete A;
  511.     if(pamiec&&b)
  512.         delete b;
  513. }
  514. void UkladPO::ZnajdzB()
  515. {
  516.     OdwrocR();
  517.     int i,j;
  518.     double pom;
  519.     for(i=0;i<No0_K;i++)
  520.     {  pom=0;
  521.         for(j=0;j<No0_K;j++)
  522.             pom+=(*RT1)[i][j]*(*b)[I_K[j]];
  523.         (*B)[i]=pom/(*d)[i+1];
  524.     }
  525.     for(i=b->Ile_Wierszy-1;i>=0;i--)
  526.     {   pom=0;
  527.         for(j=i+1;j<b->Ile_Wierszy;j++)
  528.             pom+=(*R)[i][j]*(*y)[j];
  529.         (*y)[i]=(*B)[i]-pom;
  530.     }
  531. }
  532.  
  533.  
  534. int UkladPO::Rozwiaz()
  535. {  int i,j;
  536.     double pomoc;
  537.     if(Q&&B)
  538.     {  if((Q->Bledy==BEZ_BLEDU)&&(B->Bledy==BEZ_BLEDU))
  539.         {  if((R->Bledy==BEZ_BLEDU)&&(y->Bledy==BEZ_BLEDU))
  540.             {  x=new Wektor(Q->Ile_Wierszy);
  541.                 if(x)
  542.                 {  if(x->Bledy==BEZ_BLEDU)
  543.                     {  for(i=0;i<No0_K;i++)
  544.                         {  pomoc=0;
  545.                             for(j=0;j<No0_K;j++)
  546.                                 pomoc+=(*R)[i][j]*(*y)[j];
  547.                             (*B)[i]=pomoc;
  548.                         }
  549.                         for(i=0;i<Q->Ile_Wierszy;i++)
  550.                          {  pomoc=0;
  551.                              for(j=0;j<No0_K;j++)
  552.                                 pomoc+=(*Q)[i][j]*(*B)[j];
  553.                              (*x)[i]=pomoc;
  554.                         }
  555.                     }
  556.                     else
  557.                         Bledy|=BRAK_PAMIECI;
  558.                 }
  559.                 else
  560.                     Bledy|=BRAK_PAMIECI;
  561.             }
  562.             else
  563.                 Bledy|=BLEDNE_OPERACJE;
  564.         }
  565.         else
  566.             Bledy|=BLEDNE_DANE;
  567.     }
  568.     else
  569.         Bledy|=BRAK_PAMIECI;
  570.     if(Bledy==BEZ_BLEDU)
  571.     {  ERR=SprERR();
  572.         return 1;
  573.     }
  574.     else
  575.         return 0;
  576. }
  577.  
  578.  
  579.  
  580. void UkladPO::OdwrocR()
  581. {  Macierz *Pom;
  582.     double pomx;
  583.     Pom=R->Transponuj();
  584.     RT1=Pom->Odwrotna();
  585.     delete Pom;
  586. }
  587.  
  588.  
  589.  
  590.